Indice | Previo: Aplicaciones.AgujerosNegros | Siguiente: Aplicaciones.AgujerosNegros.FenomenologiaBasica

3.2.3. Métrica de Schwarzschild

Queremos encontrar la métrica en el exterior de una distribución de masa con un total de masa M, con simetría esférica.

Por razones de simetría usaremos coordenadas esféricas $(x^0,x^1,x^2,x^3)=(ct,r,\theta,\phi)$ que en este contexto se conocen como coordenadas de Schwarzschild.

Por estar en el exterior de la masa, la ecuación de campo se reduce a la ecuación en el vacío:

$$ R_{\mu\nu}=0 $$

Infinitas métricas $g_{\mu\nu}$ son compatibles con esta condiciones. ¿Qué es entonces lo que determina la métrica? Las simetrías y condiciones de frontera:

  • La métrica debe tener simetría esférica. Esta simetría implica:

    • Que una forma simple de la métrica se obtiene usando coordenadas esféricas o coordenadas de Schwarzchild.
    • Que los coeficientes métricos deberían solo depender de la coordenada radia $r$.
    • Que la métrica es simétrica, $g_{\mu\nu}=0$ si $\mu\neq\nu$.

  • Suponemos que la métrica es estática. Esto significa que no depende del tiempo (estacionaria) pero además no cambia si se invierte el tiempo $t\rightarrow -t$.

  • La métrica es asintóticamente plana. Esto implica que a gran distancia $t\rightarrow \infty$ la métrica se reduce a la métrica de Minkowski en coordenadas esféricas:

    $$ \mathrm{d}s^2=\mathrm{d}T^2-\mathrm{d}r^2-\mathrm{d}{\ell}^2 $$ donde $T=ct$:

    $$ \mathrm{d}{\ell}=r^2\mathrm{d}\theta^2+r^2\sin^2\theta\mathrm{d}\varphi^2 $$ es el elemento de línea sobre la superficie de una esfera.

La métrica más general que cumple esas condiciones es:

$$ \mathrm{d}s^2=\mathrm{e}^{2A(r)}\mathrm{d}T^2-\mathrm{e}^{2B(r)}\mathrm{d}r^2-\mathrm{d}{\ell}^2 $$

El uso de la función exponencial tiene como proposito:

  1. Garantizar que los coeficientes de la métrica son positivos y la signatura se mantiene (métrica de un espacio Lorentziano).
  2. Facilitar las derivadas.

Nótese también que la parte tangencial de la métrica no depende de $r$ para garantizar la simetría esférica.

La condición de que la métrica debe ser asintóticamente plana implica que:

$$ \lim_{r\rightarrow \infty}A(r)=1\\ \lim_{r\rightarrow \infty}B(r)=1 $$

Tensores asociados a la métrica {tensores_schw}

Para procecer a obtener información sobre las funciones $A(r)$ y $B(r)$ debemos calcular todas las cantidades tensoriales asociadas a la métrica y reemplazar en la ecuación de campo.

Las componentes de la métrica son:

$$ \begin{array}{rcl} g_{TT} & = & \mathrm{e}^{2A(r)}\\ g_{rr} & = & -\mathrm{e}^{2B(r)}\\ g_{\theta\theta} & = & -r^2\\ g_{\varphi\varphi} & = & -r^2\sin^2\theta\\ \end{array} $$

Y las de la inversa:

\begin{eqnarray} \nonumber g^{TT} & = & \mathrm{e}^{-2A(r)}\\ \nonumber g^{rr} & = & -\mathrm{e}^{-2B(r)}\\ \nonumber g^{\theta\theta} & = & -\frac{1}{r^2}\\ \nonumber g^{\varphi\varphi} & = & -\frac{1}{r^2\sin^2\theta}\\ \end{eqnarray}

El procedimiento es bastante laborioso y para hacerlo se puede utilizar una herramienta de cálculo simbólico como Wolfram Lab Cloud. En el sitio electrónico del libro se deja al estudiante las rutinas y una guía para los comandos para realizar este cálculo.

Descargue el archivo de texto plano con las instrucciones de este enlace.

Símbolos de Christoffel {christoffel_schw}

Dado que la métrica es simétrica, los únicos símbolos de Christoffel distintos de cero serán:

\begin{eqnarray} \nonumber \Gamma^{\lambda}_{\lambda\nu} & = & +\frac{1}{2} g^{\lambda\lambda} g_{\lambda\lambda,\nu} \\ \nonumber \Gamma^{\lambda}_{\mu\mu} & = & -\frac{1}{2} g^{\lambda\lambda} g_{\mu\mu,\lambda} \end{eqnarray}

Los símbolos de Christoffel diferentes de cero serán:

\begin{eqnarray} \nonumber \Gamma^{T}_{Tr} &= & A^{\prime}\\ \nonumber \Gamma^{r}_{rr} &= & B^{\prime} \\ \nonumber \Gamma^{\theta}_{\theta r} & = & \frac{1}{r} \\ \nonumber \Gamma^{\varphi}_{\varphi r} &= & \frac{1}{r} \\ \nonumber \Gamma^{\varphi}_{\varphi\theta} &=&\cot \theta\\ \nonumber \Gamma^{r}_{TT} &=&A^{\prime} \mathrm{e}^{2(A-B)} \\ \nonumber \Gamma^{r}_{\theta\theta} &=&-r \mathrm{e}^{-2 B} \\ \nonumber \Gamma^{r}_{\varphi\varphi} &=&-\mathrm{e}^{-2 B} r \sin ^{2} \theta \\ \nonumber \Gamma^{\theta}_{\varphi\varphi} &=&-\sin \theta \cos \theta \\ \end{eqnarray}

Tensor de Riemann {riemann_schw}

El tensor de Riemann es:

$$ {R^\lambda}_{\mu\nu\kappa}=\{\Gamma^\lambda_{\mu\nu,\kappa} + \Gamma_{\mu\nu}^\delta \Gamma_{\kappa\delta}^\lambda \}-\{\nu\leftrightarrow\kappa\} $$

que se puede expandir como:

$$ {R^\lambda}_{\mu\nu\kappa}=\Gamma^\lambda_{\mu\nu,\kappa}-\Gamma^\lambda_{\mu\kappa,\nu}+ \Gamma_{\mu\nu}^\delta \Gamma_{\kappa\delta}^\lambda-\Gamma_{\mu\kappa}^\delta \Gamma_{\nu\delta}^\lambda $$

Por las simetría del tensor de Riemann y el hecho que la métrica sea diagonal las únicas componentes distintas de cero son:

$$ {R^{\lambda}}_{\mu\nu\kappa},{R^{\lambda}}_{\mu\nu\lambda} $$

es decir:

$$ {R^\lambda}_{\mu\nu\lambda}=\Gamma^\lambda_{\mu\nu,\lambda}-\Gamma^\lambda_{\mu\lambda,\nu}+ \Gamma_{\mu\nu}^\delta \Gamma_{\lambda\delta}^\lambda-\Gamma_{\mu\lambda}^\delta \Gamma_{\nu\delta}^\lambda $$

En particular los términos distintos de cero serán:

$$ {R^\lambda}_{\mu\mu\lambda}=\Gamma^\lambda_{\mu\mu,\lambda}-\Gamma^\lambda_{\mu\lambda,\mu}+ \Gamma_{\mu\mu}^\delta \Gamma_{\lambda\delta}^\lambda-\Gamma_{\mu\lambda}^\delta \Gamma_{\mu\delta}^\lambda $$

Un laborioso cálculo permite encontrar las componentes diferentes de cero no repetidas:

$$ \begin{array}{l} {R^{T}}_{rrT}=A^{\prime} B^{\prime}-A^{\prime \prime}-\left(A^{\prime}\right)^{2} \\ {R^{T}}_{\theta \theta T}=-r \mathrm{e}^{-2 B} A^{\prime} \\ {R^{T}}_{\varphi \varphi T}=-r \mathrm{e}^{-2 B} A^{\prime} \sin ^{2} \theta \\ {R^{r}}_{\theta \theta r}=r \mathrm{e}^{-2 B} B^{\prime} \\ {R^{r}}_{\varphi \varphi r}=r \mathrm{e}^{-2 B} B^{\prime} \sin ^{2} \theta \\ {R^{\theta}}_{\varphi \varphi \theta }=\left(1-\mathrm{e}^{-2 B}\right) \sin ^{2} \theta \end{array} $$

Tensor de Ricci {ricci_schw}

El tensor de Ricci se define como:

$$ R_{\mu\nu}=\sum_\lambda {R^{\lambda}}_{\mu\nu\lambda} $$

Dado que los únicos términos distintos de cero del tensor de Riemann son aquellos para los cuales $\mu=\nu$, descubrimos que el tensor de Ricci es diagonal:

$$ R_{\mu\mu}=\sum_\lambda {R^{\lambda}}_{\mu\mu\lambda} $$

Para el cálculo por ejemplo de $R_{TT}$, por ejemplo:

$$ R_{TT}={R^r}_{TTr}+{R^\theta}_{TT\theta}+{R^\varphi}_{TT\varphi} $$

Para obtener el término ${R^r}_{TTr}$ se debe usar la propiedad:

$$ {R^r}_{TTr}=g^{rr}R_{rTTr}=g^{rr}R_{TrrT}=g^{rr}g_{TT}{R^T}_{rrT} $$

Después de un álgebra laboriosa obtenemos:

\begin{eqnarray} \nonumber R_{TT} & = & -\mathrm{e}^{2(A-B)}\left[A^{\prime \prime}+\left(A^{\prime}\right)^{2}-A^{\prime} B^{\prime}+\frac{2 A^{\prime}}{r}\right] \\ \nonumber R_{rr} & = & A^{\prime \prime}+\left(A^{\prime}\right)^{2}-A^{\prime} B^{\prime}-\frac{2 B^{\prime}}{r} \\ \nonumber R_{\theta\theta} & = & \mathrm{e}^{-2 B}\left[1+r\left(A^{\prime}-B^{\prime}\right)\right]-1 \\ \nonumber R_{\varphi\varphi} & = & \sin ^{2} \theta\left\{\mathrm{e}^{2 B}\left[1+r\left(A^{\prime}-B^{\prime}\right)\right]-1\right\} \end{eqnarray}

Escalar de curvatura {escalar_curvatura_schw}

El escalar de curvcatura se calcula usando:

$$ R=g^{\mu\nu}R_{\mu\nu}=g^{00} R_{00}+g^{11} R_{11}+g^{22} R_{22}+g^{33} R_{33} $$

Al reemplazar queda:

$$ R=-2 \mathrm{e}^{-2 B}\left[A^{\prime \prime}+\left(A^{\prime}\right)^{2}-A^{\prime} B^{\prime}+\frac{2}{r}\left(A^{\prime}-B^{\prime}\right)+\frac{1}{r^{2}}\right]+\frac{2}{r^{2}} $$

Tensor de Einstein {tensor_einstein_schw}

Por último el tensor de Einstein es la traza inversa del tensor de Ricci:

$$ G_{\mu\nu}=R_{\mu\nu}-\frac{1}{2}g_{\mu\nu}R $$

Dado que ambos $g_{\mu\nu}$ y $R_{\mu\nu}$ son simétricos el tensor de Einstein también lo es:

$$ G_{\mu\mu}=R_{\mu\mu}-\frac{1}{2}g_{\mu\mu}R $$

El álgebra en este caso produce:

\begin{eqnarray} \nonumber G_{TT} & = &-\frac{2 \mathrm{e}^{2(A-B)}}{r} B^{\prime}+\frac{\mathrm{e}^{2(A-B)}}{r^{2}}-\frac{\mathrm{e}^{2 A}}{r^{2}} \\ \nonumber G_{rr} & = & -\frac{2 A^{\prime}}{r}+\frac{\mathrm{e}^{2 B}}{r^{2}}-\frac{1}{r^{2}} \\ \nonumber G_{\theta\theta} & = &-r^{2} \mathrm{e}^{-2 B}\left[A^{\prime \prime}+\left(A^{\prime}\right)^{2}+\frac{A^{\prime}-B^{\prime}}{r}-A^{\prime} B^{\prime}\right] \\ \nonumber G_{\varphi\varphi} & = &-r^{2} \mathrm{e}^{-2 B} \sin ^{2} \theta\left[A^{\prime \prime}+\left(A^{\prime}\right)^{2}+\frac{A^{\prime}-B^{\prime}}{r}-A^{\prime} B^{\prime}\right] \end{eqnarray}

Ecuación de campo y solución {ecuacion_campo_schw}

Con todos los elementos a la mano, podemos ahora preguntarnos por cuál es la condición que debe cumplir $A(r)$ y $B(r)$ para que se cumpla la ecuación de campo:

$$ G_{\mu\nu}=\frac{8\pi G}{c^4}T_{\mu\nu} $$

o bien:

$$ R_{\mu\nu}=\frac{8\pi G}{c^4}\left(T_{\mu\nu}-\frac{1}{2}g_{\mu\nu}T\right) $$

Como estamos en un punto por fuera del cuerpo, eso implica que ambos:

$$ G_{\mu\nu}=0 $$

y

$$ R_{\mu\nu}=0 $$

Adicionalmente esta última condición implica también que:

$$ R=0 $$

Tenemos un total de 9 condiciones posibles:

\begin{eqnarray} \nonumber -\mathrm{e}^{2(A-B)}\left[A^{\prime \prime}+\left(A^{\prime}\right)^{2}-A^{\prime} B^{\prime}+\frac{2 A^{\prime}}{r}\right] & = & 0\;(R_{TT})\\ \nonumber A^{\prime \prime}+\left(A^{\prime}\right)^{2}-A^{\prime} B^{\prime}-\frac{2 B^{\prime}}{r} & = & 0\;(R_{rr}) \\ \nonumber \mathrm{e}^{-2 B}\left[1+r\left(A^{\prime}-B^{\prime}\right)\right]-1 & = & 0\;(R_{\theta\theta}) \\ \nonumber \sin ^{2} \theta\left\{\mathrm{e}^{2 B}\left[1+r\left(A^{\prime}-B^{\prime}\right)\right]-1\right\} & = & 0\;(R_{\varphi\varphi})\\ \nonumber -2 \mathrm{e}^{-2 B}\left[A^{\prime \prime}+\left(A^{\prime}\right)^{2}-A^{\prime} B^{\prime}+\frac{2}{r}\left(A^{\prime}-B^{\prime}\right)+\frac{1}{r^{2}}\right]+\frac{2}{r^{2}} & = & 0\;(R)\\ \nonumber -\frac{2 \mathrm{e}^{2(A-B)}}{r} B^{\prime}+\frac{\mathrm{e}^{2(A-B)}}{r^{2}}-\frac{\mathrm{e}^{2 A}}{r^{2}} & = & 0\;(G_{TT})\\ \nonumber -\frac{2 A^{\prime}}{r}+\frac{\mathrm{e}^{2 B}}{r^{2}}-\frac{1}{r^{2}} & = & 0\;(G_{rr}) \\ \nonumber -r^{2} \mathrm{e}^{-2 B}\left[A^{\prime \prime}+\left(A^{\prime}\right)^{2}+\frac{A^{\prime}-B^{\prime}}{r}-A^{\prime} B^{\prime}\right] & = & 0\;(G_{\theta\theta}) \\ \nonumber -r^{2} \mathrm{e}^{-2 B} \sin ^{2} \theta\left[A^{\prime \prime}+\left(A^{\prime}\right)^{2}+\frac{A^{\prime}-B^{\prime}}{r}-A^{\prime} B^{\prime}\right] & = & 0\;(G_{\varphi\varphi}) \end{eqnarray}

Hay muchas posibilidades pero una de las más sencillas es combinar $\mathrm{e}^{-2A}G_{TT}+\mathrm{e}^{-2B}G_{rr}$ para producir:

$$ \frac{2 \mathrm{e}^{-2 B}}{r}\left(A^{\prime}+B^{\prime}\right)=0 $$

que implica obviamente:

$$ A'=-B' $$

Está última ecuación se puede integrar para producir:

$$ B(r)=-A(r)+C $$

¿Qué papel tiene la constante C? Lo único que hace $C$ es hacer el coeficiente de la componente $r$ de la métrica $\exp[2B(r)]=\exp[-2A(r)+2C]=\exp[-2A(r)]\exp(2C)$ un poco mayor o un poco menor en una constante, sin afectar en nada su dependencia de $r$. Ese factor adicional se puede eliminar con un cambio en las unidades de $r$, tal que en esas unidades $\exp(2C)=1$ y por tanto $C=0$. La relación resultante es entonces:

$$ B(r)=-A(r) $$

De modo que la métrica será:

$$ \mathrm{d}s^2=\mathrm{e}^{2A(r)}\mathrm{d}T^2-\mathrm{e}^{-2A(r)}\mathrm{d}r^2-\mathrm{d}{\ell}^2 $$

Con esta identidad la ecuación para $G_{TT}$ se convierte en:

\begin{eqnarray} \nonumber -\frac{2 \mathrm{e}^{-4B}}{r} B^{\prime}+\frac{\mathrm{e}^{-4B}}{r^{2}}-\frac{\mathrm{e}^{-2 B}}{r^{2}} & = & 0 \\ \nonumber -2 \mathrm{e}^{-2B}r B^{\prime}+\mathrm{e}^{-2B}-1 & = & 0 \\ \nonumber \frac{\mathrm{d}}{\mathrm{d}r}[r(\mathrm{e}^{-2B}-1)] & = & 0 \\ \end{eqnarray}

Esta ecuación se integra trivialmente como:

$$ \mathrm{e}^{-2B}=1+\frac{K}{r} $$

donde $K$ es una constante.

De aquí y por $A(r)=-B(r)$ obtenemos:

$$ \mathrm{e}^{2A}=1+\frac{K}{r} $$

De modo que arribamos a la métrica más general compatible con las condiciones del problema:

$$ \mathrm{d}s^2=\left(1+\frac{K}{r}\right)\mathrm{d}T^2-\left(1+\frac{K}{r}\right)^{-1}\mathrm{d}r^2-\mathrm{d}{\ell}^2 $$

Es posible probar con los programas provistos con la versión electrónica del libro que está métrica tiene tensor de Ricci, tensor de Einstein y escalar de curvatura cero como era de esperarse. Adicionalmente vemos que la métrica es asintóticamente plana, cumpliendo con todas las condiciones impuestas inicialmente.

Ahora bien: ¿cuánto vale la constante $K$?. Para esto nos podemos valer del principio de consistencia.

Sabemos que esta métrica en el límite de campo débil debe ser igual a la métrica e Newton:

$$ \mathrm{d}s^2=\left(1-\frac{R_S}{r}\right)\;\mathrm{d}T^2-\mathrm{d}r^2-\mathrm{d}{\ell}^2 $$

donde $R_S=2GM/c^2$ es la escala de longitud asociada a la masa $M$.

Una comparación entre estas dos métricas nos muestra que:

$$ K=-R_S $$

y la solución queda finalmente:

$$ \mathrm{d}s^2=\left(1-\frac{R_S}{r}\right)\mathrm{d}T^2-\left(1-\frac{R_S}{r}\right)^{-1}\mathrm{d}r^2-\mathrm{d}{\ell}^2 $$

3.2.4. Cantidades asociadas a la métrica

Con la forma final de la métrica obtenida podemos escribir finalmente los símbolos de Christoffel:

\begin{eqnarray} \nonumber \Gamma^{T}_{Tr} &=&\frac{R_S}{2r^2}\frac{1}{1-R_S/r}\\ \nonumber \Gamma^{r}_{rr} &=&-\frac{R_S}{2r^2}\frac{1}{1-R_S/r}\\ \nonumber \Gamma^{\theta}_{\theta r} &=&\frac{1}{r} \\ \nonumber \Gamma^{\varphi}_{\varphi r} &=&\frac{1}{r} \\ \nonumber \Gamma^{\varphi}_{\varphi\theta} &=&\cot \theta\\ \nonumber \Gamma^{r}_{TT} &=&\frac{R_S}{2r^2}\left(1-\frac{R_S}{r}\right)\\ \nonumber \Gamma^{r}_{\theta\theta} &=&-r\left(1-\frac{R_S}{r}\right)\\ \nonumber \Gamma^{r}_{\varphi\varphi} &=&-r\left(1-\frac{R_S}{r}\right)\sin ^{2} \theta \\ \nonumber \Gamma^{\theta}_{\varphi\varphi} &=&-\sin \theta \cos \theta \\ \end{eqnarray}

Las componentes del tensor de Riemann son:

\begin{eqnarray} \nonumber {R^{T}}_{rrT} & = & \frac{R_S}{r^3}\frac{1}{1-R_S/r}\\ \nonumber {R^{T}}_{\theta \theta T} = {R^{r}}_{\theta \theta r} & = & -\frac{R_S}{2r} \\ \nonumber {R^{T}}_{\varphi \varphi T} = {R^{r}}_{\varphi \varphi r} = -{R^{\theta}}_{\varphi \varphi \theta } & = & -\frac{R_S}{2r}\sin ^{2} \theta \\ \end{eqnarray}

Las ecuaciones geodésicas serán:

\begin{eqnarray} \nonumber \dot{U}_T & = & -\frac{R_S}{r^2}\frac{1}{1-R_S/r}U_T U_r\\ \nonumber \dot{U}_r & = & -\frac{R_S}{2r^2}\left(1-\frac{R_S}{r}\right)U_T^2+\\ \nonumber & & +\frac{R_S}{2r^2}\frac{1}{1-R_S/r}U_r^2 +\\ \nonumber & & +\left(1-\frac{R_S}{r}\right) r U_\theta^2 +\\ \nonumber & & +\left(1-\frac{R_S}{r}\right) r \sin^2\theta U_\varphi^2\\ \nonumber \dot{U}_\theta & = & -\frac{2}{r}U_r U_\theta+\cos\theta\sin\theta\;U_\varphi^2\\ \nonumber \dot{U}_\varphi & = & -\frac{2}{r}U_r U_\varphi-2\cot\theta\;U_\theta U_\varphi\\ \end{eqnarray}

donde $U_\mu\equiv\mathrm{d}x^\mu/\mathrm{d}\sigma$.

3.2.5. Integración numérica

Para conocer las propiedades de la métrica, podemos programar las ecuaciones de la geodésica e integrarlas numéricamente.

Lo primero es escribir las ecuaciones diferenciales:

In [1]:
def geodesica_schwarzschild(Y,s,Rs):
    #Coordenadas y sus velocidades
    T,r,q,f,dTds,drds,dqds,dfds=Y
    #Factores
    A=1-Rs/r
    v=Rs/(2*r**2)
    #ecuaciones
    from numpy import zeros_like,sin,cos
    dYds=zeros_like(Y)
    #Ecuaciones de velocidades
    dYds[:4]=Y[4:]
    #Aceleraciones
    dYds[4]=-2*v/A*dTds*drds
    dYds[5]=-v*A*(dTds)**2+\
             v/A*(drds)**2+\
             A*r*(dqds)**2+\
             A*r*sin(q)**2*(dfds)**2
    dYds[6]=-(2/r)*drds*dqds+cos(q)*sin(q)*(dfds)**2
    dYds[7]=-(2/r)*drds*dfds-2*cos(q)/sin(q)*(dqds)*(dfds)
    return dYds

Definamos las propiedades

In [2]:
# Constantes
c=3e8 # m/s
G=6.67e-11 # m^3 kg^-1 s^-2
M=1.98e30 # kg

# Radio de Schwarzschild
Rs=2*G*M/c**2

# Velocidad característica
from numpy import pi,sqrt
r0=10*Rs
v0=sqrt(G*M/r0)
df0ds=(v0/r0/c) # r df/d(ct) = v/c
P=2*pi*c*r0/v0 # Periodo típico

R_S = 2934.8 m
v/c = 0.22360679774997896 c
d(fi)/ds = 7.619149439484087e-06
P = 824657.0509063328 m

Las condiciones iniciales:

In [4]:
from numpy import pi,array
Y0s=array([
    0,r0,pi/2,0.0,
    1.0,0.0,0.0,0.9*df0ds
])

from numpy import linspace
ss=linspace(0,3*P,300)

La solución es:

In [5]:
from scipy.integrate import odeint
Ys=odeint(geodesica_schwarzschild,Y0s,ss,args=(Rs,))

Ts=Ys[:,0]
rs=Ys[:,1]
qs=Ys[:,2]
fs=Ys[:,3]

En coordenadas cartesianas:

In [6]:
from numpy import sin,cos
xs=rs*sin(qs)*cos(fs)
ys=rs*sin(qs)*sin(fs)
zs=rs*cos(qs)

Un gráfico de la solución

In [7]:
%matplotlib inline

In /Users/jzuluaga/anaconda/lib/python3.6/site-packages/matplotlib/mpl-data/stylelib/_classic_test.mplstyle: 
The text.latex.preview rcparam was deprecated in Matplotlib 3.3 and will be removed two minor releases later.
In /Users/jzuluaga/anaconda/lib/python3.6/site-packages/matplotlib/mpl-data/stylelib/_classic_test.mplstyle: 
The mathtext.fallback_to_cm rcparam was deprecated in Matplotlib 3.3 and will be removed two minor releases later.
In /Users/jzuluaga/anaconda/lib/python3.6/site-packages/matplotlib/mpl-data/stylelib/_classic_test.mplstyle: Support for setting the 'mathtext.fallback_to_cm' rcParam is deprecated since 3.3 and will be removed two minor releases later; use 'mathtext.fallback : 'cm' instead.
In /Users/jzuluaga/anaconda/lib/python3.6/site-packages/matplotlib/mpl-data/stylelib/_classic_test.mplstyle: 
The validate_bool_maybe_none function was deprecated in Matplotlib 3.3 and will be removed two minor releases later.
In /Users/jzuluaga/anaconda/lib/python3.6/site-packages/matplotlib/mpl-data/stylelib/_classic_test.mplstyle: 
The savefig.jpeg_quality rcparam was deprecated in Matplotlib 3.3 and will be removed two minor releases later.
In /Users/jzuluaga/anaconda/lib/python3.6/site-packages/matplotlib/mpl-data/stylelib/_classic_test.mplstyle: 
The keymap.all_axes rcparam was deprecated in Matplotlib 3.3 and will be removed two minor releases later.
In /Users/jzuluaga/anaconda/lib/python3.6/site-packages/matplotlib/mpl-data/stylelib/_classic_test.mplstyle: 
The animation.avconv_path rcparam was deprecated in Matplotlib 3.3 and will be removed two minor releases later.
In /Users/jzuluaga/anaconda/lib/python3.6/site-packages/matplotlib/mpl-data/stylelib/_classic_test.mplstyle: 
The animation.avconv_args rcparam was deprecated in Matplotlib 3.3 and will be removed two minor releases later.
In [8]:
import matplotlib.pyplot as plt
plt.close("all")
from matplotlib.patches import Circle
fig=plt.figure(figsize=(5,5))
ax=fig.gca()

ax.plot(xs/Rs,ys/Rs,'k-')
cRs=Circle((0,0),1,fill=None,ls='--')
ax.add_patch(cRs)

#Decoración
rang=1.2*max(abs(xs).max(),abs(ys).max())/Rs
ax.set_xlim((-rang,rang))
ax.set_ylim((-rang,rang))
ax.set_xlabel("$x/R_S$")
ax.set_ylabel("$y/R_S$")
ax.xaxis.set_label_coords(0.5,0.05)
ax.yaxis.set_label_coords(0.05,0.5)
ax.grid()
fig.tight_layout()

Figura 3.60.

Para ver una versión animada de esta gráfica vea la versión electrónica de este libro.

Podemos reunir los elementos anteriores usando el siguiente código:

In [9]:
#Condiciones iniciales
r0=20*Rs

#Solución a la ecuación
from numpy import pi,sqrt
v0=sqrt(G*M/r0)
df0ds=(v0/r0/c) # r df/d(ct) = v/c
P=2*pi*c*r0/v0 # Periodo típico
from numpy import pi,array
Y0s=array([
    0,r0,pi/2,0.0,
    1.0,0.0,0.0,0.8*df0ds
])
from numpy import linspace
Nt=300
ss=linspace(0,3*P,Nt)

#Solución
from scipy.integrate import odeint
Ys=odeint(geodesica_schwarzschild,Y0s,ss,args=(Rs,))
Ts=Ys[:,0]
rs=Ys[:,1]
qs=Ys[:,2]
fs=Ys[:,3]
from numpy import sin,cos
xs=rs*sin(qs)*cos(fs)
ys=rs*sin(qs)*sin(fs)
zs=rs*cos(qs)

#Gráfico
import matplotlib.pyplot as plt
plt.close("all")
from matplotlib.patches import Circle
fig=plt.figure(figsize=(5,5))
ax=fig.gca()

punto,=ax.plot([],[],'bo')
trayectoria,=ax.plot([],[],'k-')
cRs=Circle((0,0),1,fill=None,ls='--')
ax.add_patch(cRs)

#Decoración
rang=1.2*r0/Rs
ax.set_xlim((-rang,rang))
ax.set_ylim((-rang,rang))
ax.set_xlabel("$x/R_S$")
ax.set_ylabel("$y/R_S$")
ax.xaxis.set_label_coords(0.5,0.05)
ax.yaxis.set_label_coords(0.05,0.5)
ax.grid()
fig.tight_layout()

def animacion(i):
    punto.set_data(xs[i]/Rs,ys[i]/Rs)
    trayectoria.set_data(xs[max(0,i-100):i]/Rs,ys[max(0,i-100):i]/Rs)
    return punto,trayectoria

#animacion(0)
from matplotlib import animation
anim=animation.FuncAnimation(fig,animacion,frames=Nt,interval=50,blit=True,repeat=False)

Figura 3.61.

Una versión más controlable de la animación se provee a continuación:

In [10]:
from IPython.display import HTML
from matplotlib import rcParams 
rcParams['animation.embed_limit']=2**128
HTML(anim.to_jshtml())
Out[10]:

Indice | Previo: Aplicaciones.AgujerosNegros | Siguiente: Aplicaciones.AgujerosNegros.FenomenologiaBasica